arXiv:1501.07782vl [physics.optics] 30 Jan 2015 
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nonlinearity 

Vitaly Lutsky and Boris A. Malomed 
Department of Physical Electronics, 

School of Electrical Engineering, 

Tel Aviv University, Tel Aviv 69978, Israel 

We introduce a model of one- and two-dimensional (ID and 2D) optical media with the x^ 
nonlinearity whose local strength is subject to cusp-shaped spatial modulation, x^ ~ , with 

a > 0, which can be induced by spatially nonuniform poling. Using analytical and numerical 
methods, we demonstrate that this setting supports ID and 2D fundamental solitons, at a < 1 and 
a < 2, respectively. The ID solitons have a small instability region, while the 2D solitons have a 
stability region at a < 0.5 and are unstable at a > 0.5. 2D solitary vortices are found too. They are 
unstable, splitting into a set of fragments, which eventually merge into a single fundamental soliton 
pinned to the cusp. Spontaneous symmetry breaking of solitons is studied in the ID system with a 
symmetric pair of the cusp-modulation peaks. 
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I. INTRODUCTION AND THE MODEL 

It is commonly known that external potentials greatly 
expand the variety of stable localized objects (solitons 
and solitary vortices) which may be created by self- 
focusing nonlinearities pQ. Recently, much attention has 
also been attracted to the creation of effective nonlinear 
potentials (alias “pseudopotentials”, as they are called 
in solid-state physics 0) by means of spatial modula¬ 
tion of local strength of the nonlinearity [3]. Most works 
on this topic have been dealing with the cubic, alias x^, 
nonlinearity. In optical media, spatially nonuniform Kerr 
nonlinearity can be induced by an accordingly designed 
inhomogeneous density of nonlinearity-inducing dopants 
[5, by an inhomogeneous distribution of detuning in a 
uniform resonant-dopant density [3], or in composite me¬ 
dia assembled of different materials |5j. Similar nonlin¬ 
earity landscapes are relevant in models of Bose-Einstein 
condensates, where virtually any landscape, controlled 
by the optical Feshbach resonance [6], can be “painted” 
in space by a rapidly moving laser beam [7]. In the same 
context, the nonlinearity can be patterned by means of 
the magnetic Feshbach resonance, using magnetic lat¬ 
tices |8j. The limit case of the spatially inhomogeneous 
X® self-interaction in the one-dimensional (ID) geome¬ 
try corresponds to a very narrow strip with strong non¬ 
linearity, which is embedded into a linear host medium. 
In that case, the strongly concentrated nonlinearity may 
be described by the delta-function 0 E3 Eg. Another 
variety of the strongly localized self-focusing x^ non¬ 
linearity is a discrete linear lattice with one m or two 
nonlinear m sites embedded into it, or externally cou¬ 
pled to the lattice m- 

A more realistic model of the singular modulation of 
the cubic nonlinearity in the ID setting was introduced 
in Ref. [13], with the local strength featuring a cusp, 
X (3) ~ \x \~ a , around the singular point, x = 0. This 
model makes it possible to extend the study of the onset 


of collapse in nonlinear wave systems, starting from the 
theory developed for the uniform space M- It was also 
found, by means of analytical and numerical methods, 
that the nonlinear Schrodinger equation for wave field 
u(x, z), with the corresponding self-attractive cubic term 
\x\~ a \u\ 2 u, gives rise to stable ID solitons in the range 
of 0 < a < 1. In the limit of a = 1, the solitons dis¬ 
appear, as their amplitude vanishes, and solitons do not 
exists at a > 1, when the singularity of the nonlinear¬ 
ity modulation is too strong. Furthermore, in Ref. |T3] 
it was demonstrated that the same singular modulation 
with a < 1 allows one to emulate the action of attractive 
nonlinearities in sub-ID spaces, with effective dimension 
D = 2(1 -a)/ (2 -a) < 1. 

As mention in Ref. m too, a natural extension of 
the analysis should be the consideration of the singu¬ 
lar modulation of the quadratic (x^) nonlinearity in a 
second-harmonic-generating medium. In particular, an 
essential advantage of the consideration of the x^ non¬ 
linearity is the possibility to extend it to the 2D space, 
where any self-focusing cubic nonlinearity would imme¬ 
diately give rise to the collapse. Another advantage of 
the consideration of x® media is the fact that the well- 
elaborated poling technique m-m makes it possible to 
realize various spatially modulated profiles of the local 
strength of the x^ interactions in ID and 2D geome¬ 
tries. The technique of quasi-phase-matching [18] may 
also help to achieve this purpose, by creating spatial pat¬ 
terns of the phase mismatch m • Such patterns do not 
directly affect the local x® strength, but they determine 
the effective local nonlinearity in terms of the cascading 
limit lain], see Eq. (Jl2| below. 

Self-trapped solitary modes, pinned to a spot carrying 
strongly localized x^ nonlinearity, which may be ap¬ 
proximated by a delta-function, X^ 2 \ x ) ~ S(x) [20], and 
modes pinned to two such spots EB, embedded into the 
linear medium, were studied previously. Another variety 
of systems with the strongly localized x^ nonlinearity is 
represented by a linear lattice with one or two sites car- 
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rying the nonlinearity [22]. It is also relevant to mention 
a recent analysis of a dissipative ID system with uniform 
nonlinearity and a localized gain region (“hot spot” 
[23]). in which stable dissipative solitons are pinned to 
the “hot spot” [24] (recent reviews of dissipative solitons 
supported by locally applied gain are presented in Refs. 
[2S] and [26]). 


In this work, we concentrate on ID and 2D systems 
with singular spatial modulation of the nonlinear¬ 
ity, which is represented by coefficients \x\~ a and r -a , 
respectively. The model is introduced in Section II. In 
Section III, we report analytical and numerical results 
for the ID version of the system. It is shown that soli¬ 
tons exist at a < 1, and they vanish in the limit of 
a = 1. In fact, this existence region is broader than 
its counterpart in the model, related to the present 
one by the cascading limit, which is a < 1/2, see Eq. 
(12) below. An essential result is the stability map for 
the solitons, which contains a small instability region 
at negative values of the -mismatch coefficient. In 
Section IV the ID model is extended to include a sym¬ 
metric pair of singular-modulation peaks. In that case, 
the main result is a symmetry-breaking bifurcation m, 
which transforms spatially symmetric solitons into asym¬ 
metric ones. The 2D system is considered in Section V, 
where it is demonstrated that fundamental 2D solitons 
exist at a < 2 (although they have tiny amplitudes for 
1 < a < 2), and are stable at a < 0.5, according to 
the stability map shown below in Fig. [13] On the other 
hand, 2D vortex solitons, which are also constructed in 
Section V, are found to be completely unstable, although 
the scenario of their instability development is different 
from the one previously known for the uniform medium. 
The paper is concluded by Section VI. 


II. THE MODEL 


By means of an obvious rescaling, we fix the mismatch 
parameter at one of the three values: 

Q = 0, +1,-1. (3) 

Stationary solutions to Eqs. 0 . © with real propa¬ 
gation constant k are looked for as 

{u(x,y,z),v(x,y,z)} = {e tkz ip(x, y), e 2tkz ip(x, y)} , 

(4) 

with functions p(x,y ) and ^(x^y) (which are complex for 
vortex modes) obeying the stationary equations: 


-kip + + r a (f*'ip = 0, (5) 

-4fcV> + -Qip+ y-V = 0. (6) 

In the ID geometry, which corresponds to the planar, 
rather than bulk, waveguide with the y( 2 ) nonlinearity, 
Eqs. 0 and § amount to the following equations: 


iu z + ^ u xx + \x\ a u*v = 0, (7) 

2 iv z + ^v xx -Qv + ^\x\~ a u 2 = 0. (8) 

Accordingly, the ID version of stationary equations 0 
and (|6| is 

~k<p + kp" + \x\~ a i>p = 0, (9) 

-(4k + Q)ip+ 1^"+ = 0, (10) 


where functions p(x) and 'ip(x) are real, with the prime 
standing for d/dx. The above-mentioned cascading limit 
corresponds to neglecting in Eq. (10), which yields 


Basic models for self-guided beams in y( 2 ) media are 
well known and have been studied in detail, as summa¬ 
rized in reviews m and El- In the present work, we fo¬ 
cus on the two-wave (degenerate, alias Type-I) quadratic 
interactions, which are described by the following set of 
scaled equations for the complex fundamental-frequency 
(FF) and second-harmonic (SH) amplitudes, u and i?, in 
the presence of the spatial singular modulation of the y( 2 ) 
nonlinearity: 

iu z +-\7 2 u + r~ a u*v = 0, (1) 

2iv z + ^V 2 v - Qv + ^r~ a u 2 = 0, (2) 

where z is the propagation distance, diffraction opera¬ 
tor (1/2)V 2 acts on transverse coordinates (x,y), r = 
y/x 2 + y 2 , the asterisk (as well as symbol c.c. used below) 
stands for the complex conjugate, and real coefficient Q 
represents the SH-FF mismatch. Positive exponent a de¬ 
termines the spatial modulation of the nonlinearity co¬ 
efficient, the standard system corresponding to a = 0. 




2 ( 4 k + Q ) 




(ii) 


the substitution of which in Eq. 0 leads to the station¬ 
ary equation with the effective cubic nonlinearity: 


i I T I —2a 

- k * + 2*" + mTQ)M«= 0 - < 12 > 

An obvious corollary of Eqs. 0 , 0 and ©, © is 
that both 2D and ID exponentially localized solutions 
(solitons) may exist under conditions 


/c > 0,4/c + Q > 0. (13) 


An exception might be provided by embedded solitons , 
for which solely the former condition, k > 0, is necessary 
[28] . However, the numerical analysis has not revealed 
embedded solitons in the present model. 

Equations 0 and Q conserve the total power (alias 
Manley-Rowe invariant) US], HU), Hamiltonian, and the 
total angular momentum: 
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(14) 

(15) 

(16) 


Dynamical invariants of the ID version of the system are 
obvious counterparts of expressions (14) and (15) . 

It is relevant to note that, as mentioned in Ref. m, 
dependence P(k) for stationary solutions with wavenum¬ 
ber k and zero mismatch, Q m 0, can be found in a 
general form on the basis of scaling properties of Eqs. 
©, © and ©, {TO): 


validity of Eqs. ( fl7| ), making it necessary to produce 
P(k) relations in fully numerical form, as is done below. 


III. ONE-DIMENSIONAL FUNDAMENTAL 
SOLITONS 


A. The variational approximation (VA) 

P^ =0) {k) = fc (3/2) -«pW = °) (fc = 1), 

P r XD =(r> (k) ~ k 1 ~ a P^ ) =a " > {k = 1). (17) We start the consideration with the ID system, noting 

that Eqs. © and ([ 10 ]) for real functions <p(x) and 

However, Q / 0 breaks the scaling invariance and the can be derived from the respective Lagrangian, 


L = - 
2 


1 f +0 ° fl 

2 7-oo 12 


V ) 2 + (VO 2 


+ 


k(f 2 + (4fc + Q)^ 2 — \x\ a (p 2 ^ 


| dx, 


(18) 


Our first aim is to look for fundamental solitons with 
the help of the variational approximation (VA), adopting 
the Gaussian ansatz , with amplitudes A, B and inverse 
squared widths p, 7 of the FF and SH components m : 

<p(x) = A exp (—px 2 ) , xp(x) = Hexp (—yx 2 ) , (19) 


with total power 

P = vV ( 2 p)A 2 + 4-y/7r/ ( 2 7 )B 2 . (20) 

The substitution of the ansatz into Lagrangian ( |l 8 | ) 
yields 




[H 2 ( 8 /c + 2 Q + y)y / p V A 2 ^P/(2k V p)] — 4r 




( 21 ) 


where T is the Gamma-function. Two variational equa¬ 
tions, dL/dB = dL/d (A 2 ) = 0, which follow from La¬ 
grangian ©, produce the following expressions for the 
FF and SH amplitudes: 


_ /~7T(2fc+ p)(y + 2 p)( x a) / 2 

“ V 2 p 2r((l-a)/2) 


A 2 _ ^(Sk + 2 Q + y)( 2 k + p)(y + 2 P ) 1 Q 
4V7P[r((l-a)/2)] 2 


( 22 ) 


These expressions make sense at a < 1, predicting that 
the solitons disappear at a = 1 , as they yield A 2 (a = 
1) = B(a = 1) = 0. Below, it is shown directly that 
this happens indeed. It is relevant to mention that the 
cascading limit, which amounts to Eq. (12), admits the 
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existence of ID solitons only at a < 1/2, according to 

Ref. p~3l . 

The remaining variational equations, dL/dp = 
dL/d^f = 0 , lead to a system of coupled quadratic equa¬ 


tions for the inverse squared widths, p and 7 , where the 
amplitudes were eliminated by means of Eqs. [ 22 ] and 


p [7 + 2(2 - a)p\ - + 2 ap) = 0 , 

7 [(2a - 3) 7 - 2 p\ + 8fc([(2a - 1) 7 + 2 p] + 2 Q [(2a - 1) 7 + 2 p] = 0. 


(24) 


These equations can be readily solved numerically. 


B. The form of ID solitons around x — 0 

The singular modulation of the nonlinearity makes it 
necessary to analyze the structure of the solitons solu¬ 
tions at x 0 , similar to how it was done for the 
model in Ref. m- Accordingly, the solution is sought 
for in the form of an expansion, 

ip{x) = <P 0 - ipi\x\ * 2 ~ a + ip(x) = ipo - Ipi\x\ 2 ~ a + 

(25) 

where a < 2 is implied. It is easy to see that this ex¬ 
pansion corresponds to a local maximum of both fields 
at x = 0 (i.e., to possible soliton solutions) under the 
condition that 


<piM> > o, V’i/V’o > o. 


(26) 


A simple calculation, performed for x 0, yields the 
following results for expansion coefficients and ?/q: 


74 

2^0 

(27) 


(2 — a) (1 — a) ’ 

Vh = 

(MO 

(28) 

(2 — a) (1 — a)' 


An immediate conclusion following from Eqs. (27) and 


(28) is that condition (26) of having a maximum at x = 0 


amounts to the following inequalities: 

ipo > 0 at 0 < a < 1 , 
ipo < 0 at 1 < a < 2 . 


(29a) 

(29b) 


Actually, only in the former case, 0 < a < 1 , the solitons 
exist, while in the latter case, 1 < a < 2 , the singularity 
is too strong and cannot support solitons. This fact can 
be simply explained by the observation that the average 
value of the scaled x ^ nonlinearity coefficient in a region 
around x = 0 , \x\ < L, is determined by integral 


X 


(2)^ = J_ 


ID 


V' 


*dx 


£ —( 1 +a) 

1 — a 


(30) 


which converges at a < 1 and diverges at a > 1. Note 
that the VA equations ( 22 ) and (23) lead to exactly the 


r 


same conclusion, showing that the solitons exist solely at 
a < 1 , even if the expansion of the variational ansatz (19) 


does not exactly correspond to exact results represented 
by Eqs. ( [25] ) and ( [27] ) , ( |28| ). In principle, a more accurate 
ansatz may be taken as p(x) = A exp (—px 2_a ) ,- 0 ( x ) = 
B exp (—yx 2_a ) , to comply with Eq. (25), but the VA 
takes quite a cumbersome form in this case. 


C. Numerical results 

Using the VA-predicted wave forms as an input, it is 
straightforward to construct a family of fundamental- 
soliton solutions of Eqs. J9| and (10) by means of the 
Newton’s method. Figure [ljdisplays typical examples of 
the fundamental ID solitons, both stable and unstable. 

In Fig. [2] the family is characterized by the depen¬ 
dence of the total power, P, defined as per Eq. (14), on 


the singular-modulation exponent, <a, and on the prop¬ 
agation constant, k. The results labeled by VA are ob¬ 
tained from the Eqs. (20), ( 22 ), and (23), with p and 


7 produced by a numerical solution of Eq. (24). Figure 

[ 2 ] demonstrates that the VA provides a reasonable, al¬ 
though imperfect, accuracy, in comparison with the nu¬ 
merical findings. 

In accordance with what was said above, Fig. da) 
confirms that solitons indeed vanish at a = 1 and do 
not exist at a > 1. The same figure demonstrates that 
dependence of the family on mismatch Q [see Eq. ®>] is 
relatively weak. However, there is an essential difference 
between Q = — 1 and Q = 0 , + 1 , as concerns stability of 
the solitons, see below. 

An obvious feature observed in Fig. ^b) (for Q = +1) 
is the positive slope of the curves, i.e., dP/dk > 0, hence 
the soliton family satisfies the Vakhitov-Kolokolov (VK) 
criterion, which is a necessary condition for their stability 
m- The same conclusion pertains to the solitons found 
at Q = 0 , which is actually an exact result, according 
to Eq. On the other hand, plots P{k) for Q = 

— 1, shown in Fig. §b), exhibit a region where the VK 
criterion is not satisfied. 

To accurately check the stability of the ID fundamen¬ 
tal solitons, we have computed instability growth rates 
for eigenmodes of small perturbations added to the sta¬ 
tionary solitons (the procedure is described in Appendix 
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(a) (b) 

FIG. 1: (Color online) (a) Numerically found profiles of the FH and SH components for a stable ID soliton, with Q = — 1, 
k = 0.32, a = 0.5. (b) The same for an unstable ID soliton, with Q m —1, k = 0.28, a = 0.5. 




(a) (b) 

FIG. 2: (Color online) (a) The total power of ID fundamental solitons, P, vs. the singularity exponent, a, at a fixed propagation 
constant, k — 1, and different values of the mismatch parameter, Q. (b) P( k ) fo r Q = +1 and several fixed values of a. Both 
numerical results and their counterparts generated by the VA [see Eqs. |l9|-([24|] are displayed. 


[Ap This analysis has confirmed that the VK criterion is 
not only necessary but, as a matter of fact, also sufficient 
for the stability of the ID fundamental solitons in the 
present context. Thus, the solitons are completely stable 
for Q = +1 and Q = 0, while the boundary between the 
stable and unstable solitons for Q = —1 is shown, in the 
plane of (fc, a), in Fig. In fact, the narrow instabil¬ 

ity region at a = 0 is akin to the known narrow instability 
domain for fundamental ID solitons in the standard x^ 


system mm 

The predicted stability and instability of the ID funda¬ 
mental solitons was also verified by direct simulations. It 
has been found that, if the integral power (14) of unstable 
solitons is smaller than the power of their stable coun¬ 
terparts, the unstable solitons suffer complete decay (not 
shown here in detail). On the other hand, Fig. [^demon¬ 
strates that unstable solitons, whose integral power ex¬ 
ceeds that for coexisting stable solitons, do not decay, but 
rather transform into breathers oscillating around stable 
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(a) (b) 

FIG. 3: (a) P(k) for the fundamental ID solitons with the negative mismatch, Q = — 1, and several fixed values of the singular- 
modulation exponent, a. This plot demonstrates that there is a region with dP/dk < 0, where the VK stability criterion 
does not hold, (b) The stability map for the solitons with Q = — 1 in the plane of k and a (for Q = 0 and Q = +1, all 
the fundamental solitons are stable). The instability boundary is produced in the same form by the VK criterion and by the 
computation of stability eigenvalues for modes of small perturbations. 


solitons. 


the three fixed values ©• 


IV. THE ID MODEL WITH A PAIR OF 
SINGULAR-MODULATION PEAKS 

To introduce the model with the double-peak spatial 
modulation, we replace the underlying equations 0, © 
by more general ones, 

iu z + 1 u xx + g(x)u*v = 0, (31) 

2 iv z + kv xx ~Qv+ 1 g(x)u 2 = 0, (32) 

and select the double-peak modulation function g{pc) in 
the form of 


g(x) = \x — A\~ a + |x + A| -a , (33) 

where the separation between the peaks is 2A, while the 
mismatch parameter may be scaled, as above, to one of 


Numerical solution of the stationary version of Eqs. 


(31) and (32) reveals three types of stationary soliton so¬ 
lutions, namely, symmetric ones with maxima separated 
by distance 2A, asymmetric solitons, for which the am¬ 
plitudes are different at x = +A and —A, and twisted 
solitons, in which the FF component is antisymmetric, 
with opposite amplitudes at x = ±A, while the SH com¬ 
ponent is symmetric. Typical examples of profiles of all 
the three types of the solitons are displayed in Fig. [5] 
Similar to other models with two separated symmetric 
peaks of the local x ^ E2] and x^ [Tl] |33j [34] nonlin¬ 
earity strength , the asymmetric solitons are generated 
from the symmetric ones by the spontaneous-symmetry¬ 
breaking bifurcation m, which occurs with the increase 
of separation 2A between the peaks, if other parameters 
are fixed. The bifurcation is illustrated by Fig. [6j which 
shows the measure of the asymmetry between the local 
powers at the two peaks, 


[\u(x = A)| 2 +4Kx = A)| 2 ] - [\u(x = —A)| 2 + 4\v(x = —A)| 2 ] 
^ [| u{x = A)| 2 + 4\v(x = A) | 2 ] + [| u{x = —A)| 2 + 4\v{x = —A) | 2 ] 


(34) 


as a function of the soliton’s wavenumber, for Q = +1 (for Q = 0 and —1 the bifurcation diagrams are quite 
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(a) (b) 

FIG. 4: (Color online) A typical example of the evolution of the FF (a) and SH (b) components of an unstable fundamental 
soliton, at k — 0.28, Q — — 1 and a = 0.5, whose integral power exceeds that of a co-existing stable soliton (P > 0.04). In this 
case, the evolution leads to the formation of a breather close to the stable soliton. 
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(a) (b) (c) 

FIG. 5: (Color online) Typical profiles of the FF and SH components of stable ID solitons supported by the double-peak spatial 
modulation of the nonlinearity: (a) a twisted soliton (A = 1.17); (b) a symmetric soliton (A — 0.59); (c) a strongly asymmetric 
soliton (A = 0.31). In all the cases, Q = +1, a = 0.5, and k = 1.2. A is the respective value of the half-separation between the 
peaks. 


similar). Coordinates of the bifurcation points are sum¬ 
marized in Fig. [7| It is observed that the location of the 
bifurcation only weakly depends on mismatch Q. 

This is a supercritical pitchfork bifurcation |30], which 
destabilizes the symmetric soliton. Accordingly, the ac¬ 
curate analysis of the stability demonstrates that the 
asymmetric solitons are always stable when they exist, 
while symmetric ones are stable and unstable when they, 
respectively, do not or do coexist with the asymmet¬ 
ric solitons. Direct simulations [see Fig. [8] demon¬ 
strates that the unstable symmetric solitons sponta¬ 
neously transform into asymmetric breathers oscillating 
around coexisting stable asymmetric solitons. 

As concerns twisted solitons, they do not undergo 
breaking of the antisymmetry, and are stable for suffi¬ 


ciently large values of distance 2 A between the peaks. An 
example of the stability map for twisted solitons is pre¬ 
sented in Fig. (|9|, which demonstrates that the twisted 
solitons are stable at 2A > 1.5. In direct simulations, un¬ 
stable twisted solitons tend to spontaneously transform 
into stable symmetric or asymmetric ones coexisting with 
them (not shown here in detail). 
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FIG. 6: (Color online) Bifurcation diagrams for the spontaneous symmetry breaking of symmetric solitons supported by the 
double-peak spatial modulation of the nonlinearity, with separation 2A betwe en t he peaks and Q — +1, at fixed values of the 
singularity exponent a: (a) A = 0.43; (b) A = 0.27. Asymmetry parameter (34) is shown vs. wavenumber k of the solitons. 
Green and red lines depict stable and unstable soliton families, respectively (see the text). 




(a) 


(b) 


FIG. 7: (Color online) (a) The value of a at the bifurcation point of the spontaneous symmetry breaking of symmetric solitons, 
vs. their propagation constant k , for different values of Q and A = 0.27, as indicated in the figure, (b) The same for fixed 
k — 1.2 and different values of Q, A and a. The black line bounds the region where the numerical method does not converge 
to stationary solutions. Here, as well as in Figs. [9^b) and[~ 
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below, red lines are not shown where they overlap with blue ones, 


and the latter are not shown where they overlap with green lines. 


V. TWO-DIMENSIONAL SOLITONS AND 
VORTICES 

A. Analytical estimates 

In the 2D setting, axially symmetric solutions of sta¬ 
tionary equations <§ and (§i, with integer vorticity m 


(fundamental 2D solitons correspond to m = 0), are 
looked for, in polar coordinates (r, #), as 

<p(r, 9) = U(r)e im \ ip(r, 0) = V(r)e 2im0 , (35) 

with real amplitudes U (r) and V (r) satisfying equations 
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FIG. 8: (Color online) The evolution of peak powers of the two components of an originally unstable symmetric soliton at 
a sat 0.3, A = 0.8, Q — 1, and k = 5.5. The soliton tends to spontaneously rearrange itself into an asymmetric state close to the 
co-existing asymmetric soliton. 



FIG. 9: (Color online) The stability map for twisted solitons with different values of Q and fixed wavenumber k — 1. The 
solitons are stable at values of A exceeding those corresponding to the displayed boundaries. The region where the numerical 
method does not converge to a stationary solution is bounded by the black line (the leftmost one). 
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+ r~ a UV = 0, (36) 

-QV+\r~ a U 2 =0, 
(37) 


which can be derived from the Lagrangian [cf. its ID 
counterpart (18)]: 
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duV (dV \ 2 
dr J I dr J 


+ 


7 \ TT2 I A 7 2m _ \ 2 

k + —^ )U 2 + [4k + - 1 ^ ' T/2 


2r 2 


+ Q\ V 2 -r~ a U 2 V 


(38) 


The VA for the 2D solitons may be based on the following 
ansatz (which implies m > 0), 


U(r) = Ar m exp (— pr 2 ) , V (r) = Br 2m exp (—yr 2 ) , 


(39) 


cf. the ID ansatz (19). The substitution of the ansatz 


into Eq. (38) yields the following effective Lagrangian: 


L = j <j 4~ m B 2 m(4k + Q)')~ l ~ 2m T (2m) + -4~ m B 2 j- 2m T (2 + 2m)\ + 




>} 


-rri p—l—rn 


[fcT (1 + m) + pT (2 + m)] — 2E> (7 + 2p) 1 2m+ 2 T ^1 + 2m — ^ | . 


(40) 


where T is again the Gamma-function, cf. ID Lagrangian 
(21). While subsequent analysis can be performed in an 


explicit form for any integer vorticity m, eventually all 
the vortices with m > 1 turn out to be unstable, only 
the fundamental 2D solitons with m = 0 having a sta¬ 
bility area shown below in Fig. [l3j Therefore, explicit 
analytical results are given here for m = 0; neverthe¬ 
less, the VA predictions for vortices have been obtained 
too, and their comparison with numerical counterparts is 
presented below in Fig. [12] 

The first two variational equations 


pressions for the amplitudes, cf. similar results 
for the ID solitons: 


A = 


V(4fc + Q + 7) (7 + 2 p) 2 ~ a (k + p) 

V / 2p7r(l-a/2) 


B = 


(7 + 2P) 1 a / 2 (/c + p) 

2pT (1 — a/2) : 


Taking these results into account, the two remaining vari¬ 
ational equations, dL/dp = dL/d^f = 0, amount, simi¬ 
lar to the ID case [cf. Eq. (24)], to a pair of coupled 
quadratic equations, 


in a circle of radius R surrounding the singular point is 
[cf. Eq. pOk] 


r " 


x dr = 


-R- 


a 


(44) 


hence the singular nonlinearity-modulation profile may 
support 2D solitons at a < 2. 

It is possible to analyze the form of the solution at 
r —>> 0, adopting an expansion similar Eq. (25), which 


= 

was used in the ID case: 



[ ex- 
and 

<p(r) = + 

tp(x) = V’o-V’i?’ 2 “ + 

• (45) 


Substituting this into Eqs. ([5]) and (|6|, one can 
find 

easily 

(41) 

<£i 

2tp 0 

(46) 

^0 

To 

1 

to 


Vh = 

^0 

(47) 

(42) 

To 

1 

to 


cf. the similar ID relations (|27|) and (28). Obviously, 


Eqs. (46) and (47) corroborate that the 2D fundamental 


solitons exist at a < 2, as at a = 2 Eq. (47) dictates 


(f 2 = 0, which means that the solution degenerates into 
a trivial one with a flat SH field and zero FF component. 


| k (7 + ap) - (2 - a) p 2 = 0, 

1 ^ 7(2 - a)7 + 4 / 0(7 — on — 2 p) + Q( 7 — on — 2 p) = 0, 

(43) 

which can be solved numerically. 

As seen from Eqs. ( |4l| and (42), these results makes 
sense for a < 2 (recall in ID the result was meaningful for 
a < 1). Indeed, the average value of the x ^ coefficient 


B. Numerical results for 2D solitons and vortices 


The Newton’s method, applied to radial equations (36) 


and (37), with the VA prediction used as the initial guess, 


readily generates families of axisymmetric fundamental 
solitons. Solitary vortices can also be found in the nu¬ 
merical form. Typical examples of radial shapes of the 
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(a) (b) (c) 

FIG. 10: (Color online) Typical radial profiles of the FF and SH components of 2D solitons. (a) A stable fundamental soliton 
(a — 0.15, Q = l,k = 1); (b) an unstable fundamental soliton (a = 0.5, Q = 1, k = 1); (c) an unstable vortex soliton with 
m — 1 (a = 0.15, Q = 1, k — 1). 


fundamental and vortex (mm 1) solitons are displayed 
in Fig. [To] 

Properties of the 2D soliton families are summarized 
in Fig. 11 Note that dependence P(k) for Q = 0 can 
be found in the exact form, according to Eq. ©. it 
is worthy to note that Fig. [Tl] demonstrates that the 
accuracy of the VA is actually better for the 2D solitons 
than for their ID counterparts, cf. Fig. [2] 

For the sake of completeness, we display similar prop¬ 
erties of the vortex-soliton family with m = 1 in Fig. [T2| 
although, as said above, all the vortices are unstable. In 
this case, the VA still provides a reasonable overall ac¬ 
curacy, although the discrepancy is larger than for the 
fundamental solitons. 

The stability of the 2D solitons was analyzed by com¬ 
putations of eigenvalues for modes of small perturbations, 
as described in Appendix [Bj This analysis yields the fol¬ 
lowing results. First, stability regions for the fundamen¬ 
tal 2D solitons are shown in Fig. [l3j These solitons are 
less stable than their ID counterparts. Indeed, it was 
demonstrated above that the ID solitons have an insta¬ 
bility region at Q = — 1 only, see Fig. [3j while the 2D soli¬ 
tons have instability regions for all values of Q = —1, 0,1. 
Moreover, it is seen in Fig. [13] that the stability area is 
slightly larger for Q = — 1 than for Q = +1. Note that, 
while the 2D fundamental solitons exist up to a = 2, as 
shown above, the stability area is limited to a < 0.5. 

Direct simulations demonstrate that the instability 
transforms unstable 2D fundamental solitons into soli¬ 
tary breathers with a small or large amplitude of the 
intrinsic oscillations, as shown in Figs. [14] and [l5j re¬ 
spectively. 

On the other hand, vortex solitons are completely un¬ 
stable against eigenmodes (B3) of azimuthal perturba¬ 
tions, with J = 1 and 2 for the vortices with m = 1, 
and J < 4 for m = 2, similar to the instability of vortex 
solitons in the uniform x^ medium, which was studied 
in detail theoretically [51] and demonstrated experimen¬ 
tally [3T. However, direct simulations of Eqs. 0 and 


(J2|, which were carried out in the Cartesian coordinates, 
as shown in Fig. 16, exhibit dynamics very different from 
that observed in the uniform medium: the instability 
splits the vortex into three fragments (two large and one 
smaller), which seem as fundamental solitons. In the 
course of the subsequent evolution, the fragments do not 
separate (as they would do in the uniform medium), but 
feature irregular rotation around the x^ singularity (in¬ 
deed, the local maximum of the x^ attracts the solitons, 
as long as they exist). Then, two of them collide and 
merge into a single soliton. Eventually, the two remain¬ 
ing solitons collide twice: the first time, they bounce from 
each other, but the second collision leads to their fusion 
into a single fundamental soliton, which stays pinned to 
the x^ singularity. The soliton keeps 65% of the total 
power, being surrounded by a conspicuous field of radia¬ 
tion “debris”, that carries the entire angular momentum 
© (we have checked that the total momentum remains 
constant in the course of the evolution). 


VI. CONCLUSION 

As a contribution to the currently developing studies 
of the dynamics of solitons in effective nonlinear poten¬ 
tials, we have introduced a model of second-harmonic¬ 
generating media with the local strength of the quadratic 
nonlinearity featuring the spatial singularity, x^ ^ 
\x\~ a in ID, and, similarly, x® 00 r ~ a in 2D. The spa¬ 
tial modulation of the nonlinearity can be implemented 
experimentally with the help of the poling technique. We 
have found, using analytical and numerical methods, that 
the robust fundamental solitons, pinned to the singular¬ 
ity points, exist at a < 1 and a < 2, respectively, in 
the ID and 2D cases (the x^ counterpart of the ID sys¬ 
tem, related to it by the cascading limit, supports soli¬ 
tons in a narrower region, a < 1/2). The ID solitons 
are stable, except for a narrow domain in the parame¬ 
ter space, while the 2D fundamental solitons are stable 
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(a) (b) 

FIG. 11: (Color online) Properties of numerically found fundamental 2D solitons, and their VA-predicted counterparts, (a) 
The integral power vs. the singular-modulation exponent, a , at fixed values of mismatch Q , and k — 1. (b) The integral power 
vs. k at several fixed values of a and Q = 1. 


II 



(a) 


(b) 


FIG. 12: (Color online) The same as in Fig. 
derived from the effective Lagrangian (|40|. 


11 but for vortices with m 


1. The variational results were produced by equations 


at a < 1/2. A noteworthy fact is that the variational 
approximation provides more accurate results for the 2D 
fundamental solitons than for their ID counterparts. In 
the 2D setting, vortex solitons have been constructed too. 
They are unstable against splitting, but in a way different 
from the known instability scenario in the uniform 
medium: typically, the unstable vortex splits into three 
fragments, which eventually merge into a single soliton 
pinned to the central singularity. In the ID system, we 
have also studied the spontaneous symmetry breaking of 


solitons pinned to a pair of singular-modulation peaks, 
which gives rise to nontrivial asymmetric pinned modes 
via the supercritical bifurcation. 

As a development of the present work, it may be inter¬ 
esting to study the symmetry breaking in the 2D system 
with two or three symmetrically placed singular peaks, 
the latter configuration being an especially interesting 
one [36] . 
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FIG. 13: (Color online) The stability map for 2D fundamental solitons. 



(a) (b) (c) 

FIG. 14: An example of a breather with a small amplitude of intrinsic oscillations, generated by the instability of a 2D 
fundamental soliton with k = 1, Q = 1 and a = 0.5. (a,b) Shapes of absolute values of the FF and SH fields at evolution stages 
corresponding to z = 4 and z — 40; (c) peak values of the FH and SH fields as functions of z. 
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FIG. 15: The same as in Fig. |14| but for a breather with a larger amplitude of intrinsic oscillations, generated by the instability 
of a fundamental soliton with k = 2, Q = 1 and a = 0.8. 
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(c) 


(d) 


FIG. 16: An example of the instability-induced splitting of the vortex soliton, with m = 1, k = 1, Q = — 1 and a = 0.1, into a 
set of three fragments, which is followed by their staged merger into two fundamental solitons, and eventually into a single one, 
pinned to the x^ singularity. The initial instability of the vortex is dominated by perturbation eigenmode (B3) with J — 1. 
Panels (a), (b), (c) and (d) display, respectively, shapes of the absolute values of the FF and SH fields at the evolution stages 
corresponding to z = 20, 85, 130 and 230 . 
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Appendix A: Eigenmodes stability analysis of ID where and e^/, are real and imaginary parts 

model of the perturbations. Substituting u p and v p into Eqs. 

0, §, we arrive at the system of linearized equations 

The stability of perturbed stationary solutions can be 
evaluated by calculation of eigenvalues for modes of small 
perturbations. The perturbed solutions of Eqs. 0, § 
are introduced as 


Up = [ 1 / 5 ( 2 :) + (e V R + ie v i)} e tkz , (Al) 

v p = [i>(x) + + ie^i)} e 2tkz , (A2) 


2 2 ^ 


ie Vz - k((p 0 + e v ) + i(c p 0xx + e Vxx ) + |a;| a (ip 0 + e v )*{tp 0 + e^) 

+ e ip) + + e ^ xx ) ~ QC0O + + 2 l x l a ^° 


0, 

0 , 


(A3) 

(A4) 


where ipo and !po are stationary solutions found by means 
of the Newton’s method, po = u(x,z ) z= 0 and ^0 = 
v(x,z) z= 0, while are complex perturbations. Fur¬ 

ther, the two complex coupled linearized equation can be 
rewritten, in the matrix form, as a system of four real 
equations: 


” C<pR " 


-0 A 0 —B~ 


’ tpR " 

€(pl 


C 0 B 0 

v 

C(pl 



0 -D 0 E 

A 

tyR 

_ I _ 

z 

_D 0 -E 0 . 


_ tyl _ 


Appendix B: Eigenmodes stability analysis of 2D 
model 

The stability analysis for the 2D model was performed 
by taking perturbed solutions as 

u = [U (r) + e u (z, r, 0)] e i ( me + kz ) (Bl) 

v = [V{r)+ e v (z,r,6)}e 2i{me+kz '> (B2) 

and looking for perturbation eigenmodes with their own 
integer vorticity, J, which is independent of m: 


with definitions 

A = k - ^Z) (2) +diag(|2; i r a ('f/>oJ), (A6) 

B = diag(|a; i | _a (^o i )), (A7) 

C = -k+^D {2) +di&g(\x i \~ a (ip 0i )), (A8) 

D = diag(i|xi|“"((/5oJ), (A9) 

E = 2jfe - tr>(2) + 1 q. (A10) 

4 2 v ’ 


Replacing D^ by the differentiation mat rix and calcu¬ 
lating eigenvalues of resulting matrix (A5), the stability 
can be examined. 


i u = [r)e l ^ Xz+je ^ + (r)e l ( x * z + JG ) ? 
■v = Cj (r)e l ( Xz+je ^ + (r)e~ l ( x * z+J9 \ 


(B3) 


where A is the respective eigenvalue (that may be com¬ 
plex), instability corresponding to Im(A) 7^ 0. Numerical 
solution of the eigenvalue problem, generated by the lin¬ 
earization of Eqs. 0.0 with respect to the small per¬ 
turbations, produces the following characteristic matrix: 


J+ A B 0 ‘ 
-A J_ 0 -B 
§ 0 E+ 0 ’ 

0 -f 0 £_ 


(B4) 


where we define 
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J+ 

J_ 

E+ 

E_ 

A 

B 




k ~2 


—2k + — 
4 

2fe — - 
4 


diag( — + D^ - diag(^)(J + m) 2 

U rf 

^"diag( — )D^ + — diag( ^ )(J — m) 2 


diag( — )D^ + — diag(i)(J + 2m) 2 

n rX 

diag( — )D^ + D^ — diag(-^-)(J — 2m) 2 


>• 


diagfo a V 0i ), 

diag(r“ a £/oJ 


Here D ii) and D l2> are the first- and second-order differ- if there is eigenvalue of G with Im(A) ^ 0. 
entiation matrices. In this case the solution is unstable 


(B5) 
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